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ABSTRACT  v 

We  consider  the  vector  initial  value  problem 
ey  = f(y,t,e),  y(0)  = y°(e)  in  the  situation  when  the  m x m 
matrix  f^(y,t,0)  is  singular  with  constant  rank  k < m and 
has  k stable  eigenvalues.  We  show  how  to  determine  the 
unique  limiting  solution  YQ  of  the  reduced  problem 
f (Yq , t,0)  = 0 and  how  to  obtain  a uniform  asymptotic  expan- 
sion of  the  solution  which  is  valid  for  small  values  of  e 
on  finite  t intervals.  A numerical  technique  is  developed 
to  calculate  the  limiting  solution  and  the  results  of  some 
examples  are  compared  with  an  existing  code  for  stiff  diff- 
erential equations. 

1.  INTRODUCTION 


% 


We  consider  the  initial  value  problem 
ey  = f (y,t,e)  , y(0,e)  = y°(e) 


(1.1) 


for  m nonlinear  differential  equations  on  a finite  interval 
0 < t i T in  the  limit  as  the  small  positive  parameter  e 
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tends  to  zero.  Familiarity  with  singular  perturbation 
problems  leads  one  to  expect  that  (under  appropriate  stabil- 
ity conditions)  the  solution  of  (1.1)  would  converge  to  a 
limiting  solution  YQ  of  the  reduced  system 

f0(Y0,t)  = f(Y0,t,O)  = 0 (1.2) 

as  e -*■  0,  at  least  away  from  an  initial  boundary  layer  region 
of  nonuniform  convergence.  For  example,  in  the  classical 
Tikhonov  problem  (cf.  Wasow  (1976)),  when  the  Jacobian  F^ 
fy(y,t,0)  has  stable  eigenvalues  for  all  y and  t (the  region 
of  stability  can  be  further  restricted),  then  (1.2)  has  a 
unique  solution  YQ(t)  which  is  the  limiting  solution  of  (1.1) 
for  t > 0.  The  solution  generally  converges  nonuniformly  at 
t = 0 since  there  is  no  reason  to  expect  that  YQ(0)  = y°(0). 
Indeed,  if  f is  infinitely  differentiable  in  y and  t and  has 
an  asymptotic  expansion  in  e then  the  solution  y(t,e)  of 
(1.1)  can  be  represented  asymptotically  in  the  form 
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y(t,£)  = Y(t,e)  + n(T,e), 


(1.3) 


throughout  0 ^ t ^ T.  The  outer  solution  Y and  the  boundary 
layer  correction  II  both  have  asymptotic  expansions  in  e,  and 
II  tends  to  zero  as  the  stretched  (or  boundary  layer)  vari- 
able 

x = t/e  (1.4) 

tends  to  infinity. 

We  wish  to  consider  (1.1)  when  matrix  f (y,t,0)  is  singu- 
lar, and  in  particular  satisfies: 

Hypothesis  (H):  fy(y,t,0)  has  constant  rank  k } 0 ^ k < m 
for  all  t in  0 £ t £ T and  all  y;  its  nonzero  eigenvalues 
have  negative  real  parts  there ; and  its  null  space  is 
spanned  by  m - k linearly  independent  eigenvectors. 

In  this  case  we  will  find  that  the  asymptotic  solution  of 
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(1.1)  still  has  the  form  (1.3)  whenever  the  reduced  system 

(1.2)  is  consistent  and  solvable,  i.e.,  whenever  (1.2)  has  at 
least  one  solution.  However,  because  fg  is  singular,  (1.2) 
no  longer  has  a unique  solution  and  additional  analysis  is 
necessary  to  determine  the  unique  limiting  solution  for 

t > 0.  We  call  such  problems  "singular  singularly-perturbed 
problems".  Two  simple  scalar  examples  illustrating  some  of 
the  possibilities  are  (i)  ey  = 1,  y"(0)  = 0 and 
(ii)  ey  « 0,  y(0)  ■ 0.  For  (i),  the  reduced  problem  1=0 
is  inconsistent,  and  while  y = x = t/e  is  a solution  of  the 
form  (1.3)  we  see  that  II  becomes  unbounded  as  x + » , For 
(ii),  the  reduced  problem  0 * 0 is  satisfied  for  all  YQ,  but 
only  the  trivial  solution  Yq  * 0 is  a limit  of  the  unique 
solution  y = 0. 

Campbell  and  Rose  (1978)  studied  constant  coefficient 
singular  singularly-perturbed  problems  of  the  form 


ey  = F(e)y 


(1.5) 


and  showed  that  the  "semistability"  condition  of  Hypothesis 
(H)  is  a necessary  and  sufficient  condition  for  a limiting 
solution  to  exist  for  all  t > 0 and  all  initial  vectors  y°. 
O'Malley  (1978)  obtained  asymptotic  solutions  of  (1.1)  in  the 
almost-linear  case  when  f(y,t,0)  = F(t)y  + g(t),  assuming 
that  the  linear  reduced  system  F(t)Yg  + g(t)  * 0 is  consis- 
tent. A preliminary  study  of  nonlinear  systems  was  reported 
in  O'Malley  and  Flaherty  (1976).  Additional  work  on  singu- 
lar singular-perturbed  problems  was  done  by  Vasil 'eva  and 
others  (cf.  Vasil 'eva  (1976)  and  the  references  contained 
therein) . 

Asymptotic  solutions  with  a different  structure  than  (1.3) 
might  result  if  initial  values  were  restricted.  For  example, 
consider  (1.5)  with 


and  suppose  that  the  initial  components  satisfy  y^  = /ffy®, 
then  we  have  the  trivial  solution  for  t > 0,  but  the  boundary 
layer  behaviour  is  determined  by  the  stretched  variable 
a = t/Sc.  More  complicated  limiting  solutions  would  occur  if 
we  allowed  "turning  points"  where  the  rank  of  fy(y,t,0) 
changes  at  particular  y and  t values.  Studies  of  these 
interesting  and  difficult  problems  are  contained  in  the  work 
of  Howes  (1978)  and  Kreiss  (1978) . The  latter  also  contains 
numerical  methods . Two  simple  scalar  examples  of  such  prob- 
lems are  ey  = -y3  + ey  and  ey  = (t  - l/2)y,  where  the  ranks 
of  fy(yQ,t,0)  change  at  y = 0 and  t = 1/2,  respectively. 

In  Section  3 of  this  paper  we  develop  asymptotic  expan- 
sions for  the  outer  solutions  Y(t,e)  of  a special  class  of 
singular  singularly-perturbed  problems  and  in  Section  4 we 
consider  more  general  problems.  Some  preliminary  linear 
algebra  is  presented  in  Section  2.  Expansions  for  the  boun- 
dary layer  correction  H(xfe)  and  a proof  of  the  asymptotic 
validity  of  our  solutions  have  also  been  obtained  and  will  be 
reported  in  O'Malley  and  Flaherty  (1978).  In  Section  5 we 
develop  a numerical  procedure  for  calculating  the  limiting 
solution  YQ(t)  which  is  based  on  the  expansion  of  Section  4 
and  in  Section  6 we  apply  this  procedure  to  some  examples 
and  discuss  the  results. 

Our  primary  motivation  for  this  work  is  the  need  to  dev- 
elop numerical  procedures  for  singularly-perturbed  (or  stiff) 
two-point  boundary  value  problems.  However,  our  results 
should  be  applicable  to  initial  value  problems  in  power  gen- 
eration and  distribution  systems,  biological  and  chemical 
reactions,  and  electrical  networks.  A new  application  is 


ill-conditioned  nonlinear  optimization  problems  (cf.  Boggs 
and  Tolle  (1977)). 


2.  ALGEBRAIC  PRELIMINARIES 


We  shall  attempt  to  find  an  asymptotic  solution  of  (1.1) 
in  the  form  given  by  (1.3).  The  decay  of  il  as  t + “ implies 
that  the  outer  solution  Y(t,e)  should  satisfy 

ey  = f(y,t,e)  (2.1) 


as  a power  series  in  e,  i.e., 

00  • 

Y(t,e)  - -I  Y.(t)eJ.  (2.2) 

j=0  J 

Under  Hypothesis  (H)  we  are  guaranteed  that 

f Oy (y » t)  = fy(y.t.O)  = ~ (y,t,0)  (2.3) 

can  be  put  into  its  reduced  echelon  form  by  an  orthogonal 
matrix  E(y,t).  Golub  (1965)  discussed  a numerical  procedure 
for  obtaining  E by  performing  a sequence  of  k Householder 
transformations.  The  differentiability  of  E follows  that  cf 
fQ  under  the  constancy  of  rank  condition  (cf.  Golub  and 
Pereyra  (1976)).  We  partition  E as 


where  Ej  is  k x m,  E2  is  (m-k)xm,  and 

E,  £„  =0 

2 Oy 

In  addition, 


Ef  E = 
Oy 


S 

0 


X 

0 


(2.5) 


(2.6) 


where 


s = E1  Vi 


x - E1 V2  * 


(2.7) 


Hypothesis  (H)  guarantees  that  S has  k stable  eigenvalues. 

We  note  that  Clasen  et  al  (1978)  used  such  constant  ortho- 
gonal matrices  E to  integrate  stiff  problems,  while  O'Malley 
(1978)  used  time  dependent  matrices  for  almost  linear  prob- 
lems . 


The  orthogonality  of  E further  implies  that 

E.E^  = 0 , E,E?  = I,  , E.E^  = I . , and  (2.8) 

12  ’ll  k 2 2 m-k 

T T 

ElV  E2E2  - . 

where  I is  the  m x m identity  matrix.  Using  the  last  rela- 
m 

tion  we  introduce  the  complementary  orthogonal  projections 


P = EjEj  , Q = E^E2 


(2.9) 


which  provide  a direct  sum  decomposition  of  m-space  with  Q 

T T 

projecting  onto  N(f  ),  the  null  space  of  f , and  P project- 
ing onto  R(f Qy) , the  range  of  fQ  . 

3.  A SPECIAL  PROBLEM:  E(y,t)  = E(t) 

In  this  section  we  examine  special  problems  (1.1)  when  the 
orthogonal  matrix  E(y,t)  introduced  in  Section  2 is  indepen- 
dent of  y.  This,  of  course,  includes  the  nearly  linear  prob- 
lems where 

f(y,t,0)  * F(t)y  + g(t) 

and  "classical"  singular  perturbation  problems  having  the 


c?i  = fi(yi»y2’t)  + esi(yi»y2*t»e) 

ey2  = eg2(y1,y2,t,e)  . 

Here,  y^  is  a k-vector,  y2  is  an  (m-k)-vector , and  Bfj/Syj  is 


(3.1) 


of  rank  k.  In  this  case  E 


m 


We  define 
z * E(t)y, 

and  further  partition  z like  E,  i.e.. 


'2  J 
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L E,y  J 


(3.2) 


Introducing  (3.2)  into  (1.1)  gives  the  following  system  for 
z: 


ezx  = h1(z1,z2,t,e)  , z^(0)  = Ei(0)y° 

(3.3a) 

z2  = h2(z1,z2,t,c)/e  , z2(0)  « E2(0)y° 

(3.3b) 

where. 

h^  = E^f(ETz,t,e)  + eE.ETz  , i = 1,2. 

(3.4) 

We  have  divided  (3.3b)  through  by  z since 

h2(Zl,z2,t,0)  = 0 

(3.5) 

necessarily  follows  if  the  reduced  system  (1.2)  for  (1.1)  is 
consistent.  This  is  because 


3h. 


(zlfz2,t,0)  = E2f0y(ETz,t)E^  = 0 


1,2 


upon  use  of  (3.4)  and  (2.5).  Thus,  h2(z1,z2,t,0)  is  a func- 
tion of  t only.  However,  the  reduced  system  (1.2)  implies 
the  corresponding  reduced  system 

hi(z1,z2,t,0)  - 0 , i = 1,2 

for  (3.3).  Hence  h^Zj.z  ,t,0)  must  be  the  trivial  function 
of  t on  0 ^ t S T for  any  and  z2,  otherwise  (1.2)  would 
have  no  solutions.  Tikhonov's  results  apply  to  (3.3)  because 
his  stability  condition  that 


■ 
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(3.7) 


Zl(t,e)  = Z j(t, e)  + AjCt.e) 
z2(t,e)  = Z2(t,e)  + eA2(x,e) 

(cf.  O'Malley  (1974)),  where  the  outer  solution  (Zj,Z2)  and 
the  boundary  layer  correction  each  have  power  series 

expansions  in  e with  the  boundary  layer  correction  decaying 
to  zero  as  T = t ■/£■*■  00 . 

Since  the  outer  solution  provides  the  asymptotic  solution 
for  t > 0,  we  must  have 

ez:  = h^Zj.Zj.t.e)  , Z2  = h2(ZltZ2,t,e)/e  (3.8) 

satisfied  as  power  series 


Z (t.o  - I Z (t)eJ  , i = 1,2, 
3-0  LJ 


(3.9) 


in  e.  The  leading  term  must  necessarily  satisfy  the  limiting 
problem 


hl^Z10’Z20,t:’0^  = 0 


(3.10a) 


Z20  * h2  (Z10»Z20>t»O)  ' z20(O)  - E2(0)y°(0).  (3.10b) 

e 

Its  unique  solution  is  obtained  since  (3.6)  and  the  implicit 

function  theorem  imply  that  the  algebraic  equation 

hj (Zj o » Z20 » = 0 can  be  uniquely  solved  for  the  k-vector 

Zl0(t)  - *(Z20(t),  t ) (3.11) 

leaving  the  nonlinear  (m-k)  th  order  initial  value  problem 
3h 

Z20  " “§T  (*(Z20’°*  Z20*t*O)  ’ Z20(O)  " E2(0>y°(0)  (3’12) 


A 


for  Z . We  shall  assume  that  the  unique  solution  to  (3.12) 


continues  to  exist  throughout  0 £ t £ T.  Note  that  the  re- 
duced system  (1.2)  implied  that  both  hj  = 0 and  h2  = 0 along 


(Zi0.Z20.t,O).  but  it  did  not  provide  equation  (3.12)  needed 


to  uniquely  determine  the  limiting  outer  solution  ^10»z2o^' 
Higher  order  terms  in  (3.9)  satisfy  linear  problems 
3h  3h 

37^  (Z10’Z20’t,O)Zlj  + TT2  (Z10’Z20,t,O)Z2j  = 8rj-l(t) 

32h  32h  (3>13) 

• 9 9 . 

Z. 


2J 


3z  3e  (Z10,Z20,t,O)Zlj  + 3z  3e 


(Z  Z ,t,0)Z  . + 
10  20  2J 


82,i-ft) 
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with  the  g.  . (t)’s  determined  by  lower  order  terms  in  the 

1»  J— 1 


outer  expansion.  One  solves  the  first  equation  for  Z ^ j as  a 


linear  function  of  z2j  * an<*  c^en  t^le  resulting  linear  differ- 


ential equation  for  z2j • Thus,  the  outer  expansion  (3.9)  can 


be  uniquely  generated  termwise  in  0 S t £ T up  to  a knowledge 
of  the  initial  value  of  the  boundary  layer  correction  compo- 


nent A2(0,e). 


The  boundary  layer  correction  is  obtained  by  noting  that 
both  (Zj,z2)  and  (Zj,Z2)  satisfy  the  differential  equations 
(3.3).  Hence,  using  (3.8)  in  (3.3)  we  have 
dA, 


= h1(Z1(er,e)  + A^x.e)  ,Z2(ex,e)  + eA2(x,e) ,ex,e) 

-h  (Z  (ex,e),Z  (ex,e) ,ex,e)  (3.14) 

dA  1 1 2 

-jjt  = [h  (Z  (ex,e)  + A1(x,e),Z2(ex,e)  + eA2 (x, e) , ex, e) 
-h2(Z1(ex,e),Z2(ex,e) ,ex,e) ] . 


We  require  Aj  and  A^  to  decay  as  x + ® and  satisfy  the  ini- 
tial condition 


Aj (0,e)  = E1  (0)y° (0)  - Z^O.e) 


(3.15) 


Taking 


'jSo  VT)‘J  1 * 1-2 


(3.16) 


we  find  that  the  leading  term  Ajg  must  satisfy  the  nonlinear 
initial  value  problem 
dA,  „ 


IT  = hi(zio(0)  + aio(t)’Z2o(0)’°’0) 

h^z^O.z^O.o.o)  , 


(3.17) 


A (0)  = E (0)y°(0)  - Z (0) 
10  1 10 


This  problem  has  a unique  exponentially  decaying  solution 

3h  ^ 

A (t)  since  (3.6)  implies  that  t — (z  ,z  ,t,0)  has  stable 

10  dzl  1 2 

eigenvalues  for  all  arguments.  Knowing  A we  can  calculate 

A _ and  successive  terms  in  (3.14).  The  details  of  this  cal- 
20 

culation  are  omitted  here  as  they  will  be  reported  elsewhere 
(cf.  O'Malley  and  Flaherty  (1978)). 

The  asymptotic  validity  of  the  expansion  (3.9)  follows 
from  Tikhonov's  theorem  (cf.  Wasow  (1976)  or  Vasil 'eva  and 
Butuzov  (1973)).  Returning  to  the  original  variables,  we 
have  found  a unique  asymptotic  solution  of  the  form  (1.3) 
with  the  outer  solution  given  by 

Y(t, e)  = ET(t)Z1(t,e)  + E^(t)Z2(t,e) 

and  with  the  exponentially  decaying  boundary  layer  correction 
given  by 

T T 

n(f,e)  = E, ( et) A (t, e)  + eE  (ex)A  (x,e). 

11  2 2 

The  result  will  even  be  valid  for  all  t 2 0 provided  that  Z 

20 

decays  exponentially  as  t + ® (cf.  Hoppensteadt  (1966)). 

4.  THE  ORIGINAL  PROBLEM 

We  now  return  to  the  original  problem  where  the  orthogonal 
matrix  E can  depend  on  y as  well  as  t.  As  noted  in  Section 
2,  the  outer  solution  (2.2)  should  satisfy  the  system  (2.1) 


ft 


as  a power  series  in  e for  t S 0.  The  leading  term  in  the 
expansion  will  satisfy  (1.2)  and,  for  j ^ 1,  fy(YQ,t,0)Yj 
will  be  successively  d termined  by  the  preceding  Y^, 
i = 0,1,...,  j-1.  This  fails  to  uniquely  determine  the  Y^'s 
since  fy(YQ,t,0)  has  rank  k < m.  We  shall  instead  find  them 
as  solutions  of  differential  equations.  To  this  end,  we 
differentiate  (2.1)  with  respect  to  t to  obtain. 


fy(Y,t ,e)t  + ft(Y,t,e)  = eY  (4.1) 

and  use  (1.2)  and  (4.1)  together. 

We  define  E^t)  = E1Cy0(t),t)  and  let  E2Q(t) ,PQ(t) , and 

QQ(t)  be  analogously  defined.  From  (2.8)  and  (2.9)  we  see 

that  P„  + Qn  = I ; thus,  we  can  write 
0 ^0  m 


V + V 


(4.2) 


and  seek  to  obtain  equations  for  P^Y  and  Q^Y.  In  particular, 
(4.1)  and  (4.2)  imply 

EI0fy«,t,e)CP0i  * Q0-i)  - E10(-ft  * A. 

Using  the  stable  matrix 


so  - E10fy«0.C-°)E10 
and  (2.9)  we  have 

V = — Aq ( [fy (Y, t, e)  - fy(Y0,t,o)]P0Y 
+ fy(Y,t,e)Q0Y  + f t (Y, t, e)  - eY> 

where, 

Ao  ■ EioVEio' 

From  (2.1)  we  have 

QqY  = Qof(Y,t,e)/e. 

Using  (4.3)  and  (4.6)  in  (4.2)  we  find 


(4.3) 


(4.4) 


(4.5) 


(4.6) 


Y = - A0Ufy(Y,t,e)  - f (Yq , t ,0) 3 t + f,  (Y,t,e) 


- eY}  + B Q f(Y,t,e)/e 


where 


B = 1 - A f . (4.8) 

0 m 0 Oy 

It  may  be  useful  to  note  that  BQ  is  a projection  with 

BP  = °»  BO  = B , and  B A„  = 0. 

0 0 0 0 0 0 0 

Setting  e = 0 in  (4.7)  yields  the  limiting  nonlinear  equa- 


Yq = -AQf  (YQ,t,0)  + B0Q0fe(Y0,t,O) 


(4.9) 


We  note  that  the  term  QQf 0y(YQ, t)Y1  is  missing  since  Q0f0y=O 
upon  use  of  (2.5)  and  (2.9). 

In  order  to  obtain  further  coefficients  Y.  it  is  necessary 

J 

to  consider  the  coefficients  of  higher  powers  of  e in  (4.7). 


Thus,  setting 

GO 

f (Y,t,e)  - Z f . (Y, t) 


(4.10) 


and  using  the  expansion  (2.2)  for  Y implies  that 
f<Y,t,e>  - £0(Y0,t)  ♦ e[fl(Y0.t)  * tjyBj.OY,] 

* s2I£2(Y0,t)  * £oy(Yo,t)Y2  + fly(Y0.t)YI 

* 5 “oyy'V^W  + °(e3) 

together  with  corresponding  expansions  for  ft(Y,t,e)  and 
fy(Y,t,e).  The  coefficient  of  e in  (4.7)  then  provides  the 
following  nonlinear  equation  for  Yj 

V-V£lt<V‘>  * f0ty(Y0’t)Yl  * !VVC) 
*foyy(Vt,Y1HVoiV.(Vc,,-Yo) 

* Vo  {£2<Yo-t)  * Wt,Yl  * 


Except  for  the  final  quadratic  term,  this  is  a nonhomogeneous 
linearization  of  the  equation  for  Yq  . Higher  order  terms 
Yj , j>2,  satisfy  linear  differential  equations  with  success- 
ively known  nonhomogeneous  terms. 

We  note  that  it  may  be  advantageous  to  obtain  differential 
equations  for  the  successive  terms  Yj  of  the  outer  expansion 
even  in  the  special  case  (see  section  3)  where  E(y,t)  is 
independent  of  y.  In  that  case,  we  solved  the  nonlinear 
algebraic  equation  (3.10a)  for  Z as  a function  of  Z20, 
followed  by  a nonlinear  initial  value  problem  (3.10b)  for 

It  may  often  be  numerically  simpler  to  solve  the  ini- 
tial value  problem  (4.9)  for  YQ(t)  and  those  for  later  terms 
successively.  We  have  not,  however,  fully  explored  both 
possibilities . 

We  will  have  to  assume,  of  course,  that  the  nonlinear 
initial  value  problems,  (4.9)  and  (4.11),  for  YQ  and  Yj  have 
solutions  on  0 5 t 5 T.  Moreover,  since  (1.2)  and  its  time 
derivative  (4.1)  are  built  into  (4.7),  consistency  with  (1.2) 
at  t = 0 implies  consistency  of  the  outer  expansion  for 
t > 0.  If  consistency  failed,  the  form  (1.3)  of  the  solution 
would  be  inappropriate.  Thus,  using  (1.2),  (2.1),  (2.2),  and 
(4.10),  we  must  have 

f (Y  (0),0)  = 0, 

foy(YQ(0),0)Y1(0)  = YQ(0)  - f 1(Yq(0) ,0)  , (4.12) 


These  equations  always  have  a solution  under  our  assumptions. 
For  example,  in  the  second  equation  we  must  have 
Yq(0)  - f 1(Yq(0) ,0)  in  the  range  of  f Qy (YQ(0) ,0) . Recall, 
however,  that  1^  = PQ  + Qq  provides  a direct  sum  decomposi- 
tion of  m space  with 


R(QJ  = N(f  *(Yn(0)f0))  and  R(Pn) 


R<f0y(Y0(O)»O))' 


‘O'  - '-0y"Os-' "'*0' 

Thus,  the  second  of  (4.12)  will  be  automatically  satisfied 

since  Q f =0  implies  that  Q [Y  (0)  - f (Y  (0),0)]  = 0. 

o oy  oo  10 

Because  f has  rank  k,  k components  of  YQ(0)  are  deter- 
mined as  a function  of  the  remaining  m-k  components. 

Indeed,  we  could  attempt  to  solve  (1.2)  for  Eio(0)Yq(0)  in 
terms  of  E20(0)YQ(0)  since  SQ(0)  (cf.  (4.3))  is  nonsingular. 
Likewise,  for  j >0,  termwise  determination  of 
f Oy (Yq (0) »0)Yj (0)  implies  that  of  EiQ(0)Yj(0)  (by  an  argument 
similar  to  the  one  preceding  (4.3)).  Thus,  EiQ(0)Yj(0),  or 
PQ(0)Yj(0),  is  determined  termwise  while  E2o(0)Yj(0),  or 
QQ(0)Yj(0),  may  be  specified.  The  purpose  of  the  boundary 
layer  correction  is  to  compensate  for  the  jump  in 
Pq(0)(Yj(0)  - y 9)  and  to  specify  the  values  of  QQ(0)Yj(0), 
j * 0. 

Once  again,  the  representation  (1.3)  and  the  fact  that  the 
differential  equation  (1.1)  is  satisfied  by  both  y and  Y 
imply  that  the  boundary  layer  correction  n(x,e)  must  satisfy 
the  nonlinear  equation 


= f(Y(er,e)  + n(x, e) , ex, e)  - f (Y(ex,e) ,ex,e) 
as  a power  series 


n(x, e)  - I n.(x)cJ 

j=0  J 

in  e and  decay  to  zero  as  x 
n(0,e)  - y°(e)  - Y(0,e). 


Moreover, 


(4.13) 


(4.14) 


(4.15) 


The  details  of  the  calculation  of  the  boundary  layer  cor- 
rection and  a proof  of  the  asymptotic  validity  of  the  solu- 
tion are  omitted  here  and  will  be  presented  in  O'Malley  and 
Flaherty  (1978).  We  summarize  our  findings,  however,  in  the 
following  theorem. 


Theorem:  Consider  the  initial  value  problem 


ey  = f (y , t,  e)  , y(0)  = y°(e) 

for  an  m-vector  y as  e -*•  0+.  Assume  that: 

(i)  f is  infinitely  differentiable  in  y and  t and  f and 

y°(e)  have  asymptotic  series  expansions  in  powers  of  e. 

(ii)  There  exists  an  infinitely  differentiable  orthogonal 
matrix  E(y,t)  for  all  y and  for  t in  the  interval 
0 2 t ^ T such  that  E(y,t)f  (y,t,0)  is  row-reduced  and 
of  rank  k,  0 ^ k < m.  Moreover,  partitioning  E after 
its  first  k rows  as  in  (2.4),  we  have 

Efy (y , t,0)ET  = 

T 

where  S = E^f^Cy, t,0)Ej  is  a stable  matrix  for  all 
values  of  y and  t . 

(iii)  The  nonlinear  system 


f (Yq (0) ,0,0)  = 0 

Q(Y  (0) ,0) [y 0 (0)-Y  (0) 

0 o 

+ f f(Y  (o)  + n (T),o,o)dT]  = o 

JO  0 

0 

can  be  uniquely  solved  for  Y^(0). 
decaying  solution  of 


(4 .16a) 

(4.16b) 

Here,  II  (t)  is  the 


dn 


= f(YQ(o)  + n0(x),o,o)  - f (Yfl(0) ,o,o) , 

n (0)  = y°(0)  - Y (0). 

0 0 


(iv)  The  matrix 


1 - coEio(0)VV0)’0) 

is  invertible  for  a particular  matrix  C . (This 
insures  that  Yj(0)  may  be  uniquely  determined.) 

(v)  The  initial  value  problems  (4.9)  and  (4.11)  have 


solutions  on  the  interval  0 ^ t ^ T. 


Then,  the  initial  value  problem  (1.1)  has  a unique  solution 
of  the  form 

y(t,e)  = Y(t,e)  + n(x,e). 

Some  comments  on  this  theorem  are  listed  below. 

(i)  Hypothesis  (ii)  is  guaranteed  by  out  earlier  Hypo- 
thesis (H) . 

(ii)  The  theorem  is  easily  obtained  from  Tikhonov's  result 
if  E(y,t)  is  independent  of  y.  It  is  considerably 
simplified  if  only  ^(Y  (0),0),  and  thereby  Q(Yq(0),0), 
is  independent  of  YQ(0).  In  this  case  (4.16b)  reduces 
■>'  to  the  linear  equation 

Q0(O)[y°(O)  - Yq(0)]  = 0 • (4.17) 

and  (4.16a)  becomes  a nonlinear  equation  for 
P0(O)Yq(O).  It  can  be  further  shown  that  the  inverti- 
bility  condition  of  Hypothesis  (iv)  then  will  be  auto- 
matically satisfied. 

(iii)  Higher  order  terms  follow  without  complication  under 
these  hypotheses. 

(iv)  Vasil 'eva  (1976)  considers  such  problems  under  a list 
of  ten  hypotheses,  generally  paralleling,  but  more 
restrictive  than  ours.  Her  most  critical  assumption 
involves  the  existence  of  a k-dimensional  manifold  of 
decaying  solutions  for  n^(x)  which  can,  more  or  less, 
be  stated  in  the  form 

E20(o)n0(t)  = *(E10(o)n0(T)) 

for  a particular  function  $ and  for  all  x and  Y (0) . 

0 

At  x = 0,  we  would  have 

E20(O)[y°(O)  - YQ (0) ] = *(E10(O)[y°(O)  - YQ(0)]) 


r 


where  YQ(0)  must  also  satisfy  the  reduced  equation  at 
t = 0.  This  analog  of  Hypothesis  (iii)  should 
uniquely  determine  Y^(0)  so  that  the  resulting  II Q (0) 
lies  on  the  manifold  of  initial  values  corresponding 
to  decaying  solutions  of  Hq(t). 

5.  NUMERICAL  ALGORITHM 

We  have  developed  an  algorithm  based  on  the  asymptotic 
analysis  of  Section  4 to  calculate  the  leading  term  Y^Ct)  in 
the  outer  solution.  For  most  problems  it  is  possible  to 
calculate  numerical  solutions  without  explicitly  identifying 
a small  parameter  e;  thus,  we  consider  initial  value  problems 
in  the  form 

y = f(y,t,e)  5 f(y,t,e)/e  , y(0)  = y‘(c).  (5.1) 

The  e,  although  shown  in  (5.1),  is  to  be  regarded  as  uniden- 
tifiable. However,  if  the  actual  limiting  solution  is 
desired,  the  evaluation  of  y in  (5.1)  causes  overflow,  or 
the  order  e terms  in  f are  close  to  the  unit  roundoff  of  the 
computer  relative  to  the  order  unity  terms  in  f,  then  a value 
of  e can  be  identified  and  the  initial  value  problem  can  be 
written  in  the  form  of  equation  (1.1).  The  actual  computer 
code  is  capable  of  handling  both  cases,  and  all  that  the  user 
need  do  is  define  f as  in  (5.1)  or  f as  in  (1.1). 

The  algorithm  consists  of  two  parts:  (i)  calculating  the 
initial  conditions  YQ(0)  for  the  outer  problem  and  (ii) 
integrating  the  differential  equation  (cf.,  (4.9))  for  YQ(t). 
We  first  describe  the  integration  procedure. 

The  differential  equation  (4.9)  for  YQ  is  not  stiff; 
hence,  any  good  code  for  integrating  non-stiff  initial  value 
problems  may  be  used.  We  use  both  the  Adams'  methods  that 
are  incorporated  into  the  Hindmarsh  (1974)  version  of  Gear's 


code  and  the  IMSL  version  of  the  Bulirsch  and  Stoer  (1966) 
extrapolation  procedure.  Both  of  these  codes  require  the 
evaluation  of  YQ  as  a function  of  YQ  and  t,  and  we  accom- 
plish this  as  follows: 

(i)  Calculate  E(YQ,t)  by  decomposing  f (Y0,t,e).  It  is  not 
necessary  to  set  e to  0 unless  e has  been  explicitly 
recognized  and  the  actual  limiting  solution  is  desired. 
Golub's  (1965)  procedure,  which  uses  a sequence  of  k 
Householder  transformations  with  column  pivoting,  is 
used  to  obtain  E.  At  the  v th  step,  1 < v < k,  of  this 
procedure  we  have 


E (Y0,t)?y(Y0,t,e) 


where,  U is  v x v and  upper  triangular,  V is 
v x (m-v),  and  K is  a permutation  matrix  due  to  the 
column  pivoting.  The  procedure  terminates,  and  the 

A 

rank  k of  f is  determined,  when  the  maximum  available 

pivot  element  in  W is  small  relative  to  the  diagonal 

(k) 

elements  of  U.  We  then  have  E **  E . The  decomposi- 
tion is  not  performed  at  every  time  step,  but,  rather 
a test  is  used  to  determine  if  E has  changed  by  too 
much.  Thus,  the  same  matrix  E may  be  used  for  several 
time  steps  or,  when  E is  constant,  for  the  entire 
integration.  If  at  any  stage  of  the  computation  the 

A 

rank  k of  f changes,  a turning  point  has  probably 
been  encountered,  and  the  integration  is  terminated, 

(ii)  Partition  E into  Ej  and  E2  as  in  (2.4).  Calculate 

Q = E^E2  and  S = Ejf  (Y0,t,e)E*. 

(iii)  Calculate  QYfl  = Qf(YQ,t,e)  and  b = -Ej [f t(Y0, t, e) 

+ fy(YQ ,t, e) (QYq) ] . When  e is  explicitly  recognized 
QYq  is  calculated  as  QfE(Yo,t,0) . 


(iv)  Solve  S(EJY0)  = b for  EjYg  by  Gaussian  elimination  and 
calculate 

■ £I(EiV- 

(v)  Calculate  Y = PY  + QY  . 

0 0 0 

We  now  turn  to  the  calculation  of  the  initial  conditions 
Yg(0)  for  the  outer  problem.  This  is  a difficult  task  when 
E2  depends  on  y.  It  requires  the  solution  of  the  nonlinear 
system  (4.16)  and  the  computation  of  the  boundary  layer  solu- 
tion II0(t),  which  itself  depends  on  YQ  (0) . It  is  possible 
that  the  integral  in  (4.16b)  may  be  adequately  approximated 
by  a very  crude  quadrature  rule,  which  would  greatly  simplify 
the  computation.  Miranker  (1973)  has  successfully  used  such 
a technique  on  stiff  problems,  but  we  have  not  as  yet 
explored  this  possibility.  Our  code  has  only  been  implemen- 
ted for  problems  where  E2  is  independent  of  y;  thus,  when  e 
is  not  explicitly  recognized  YQ(0)  is  determined  as  the  solu- 
tion of 


f (Yq (0) ,0,e)  = 0 

yy0)  - y°<e)i  = o. 


(5.2) 


A Newton  like  iteration  scheme,  which  closely  parallels  the 
computation  of  YQ(t)  is  used  to  solve  this  nonlinear  system. 
The  procedure  is  outlined  below. 

(i)  Select  an  initial  guess  for  YQ(0),  e.g.,  * y 

and  set  p = 0. 

(ii)  Calculate  by  decomposing  f (X^,0,e).  This  is 

performed  as  in  step  (i)  of  the  procedure  for 

calculating  YQ.  If  Ej  is  independent  of  y then  this 
step  need  only  be  performed  once. 

(iii)  Calculate  and  as  in  step  (ii)  of  the  previ- 

ous procedure. 


(iv)  Calculate  q(w+1)  - Q(ij)  (y°-X(u) ) and 


» 


- -E 


(y)ri/V(u) 


[ f (X  ,0, e)  + f (XVM'0,E)q 


(y), 


.(y+i) 


y 

(v)  Solve  (X(y+l)-  x(u))] 

-X^)  and  calculate  = 

-x(li))] 


for  E^(X^  + 1>. 

(yKTri,(y)  ,„(y+i). 


(e;m/)‘[e 


(xv 


(vi)  Set  X(y+1)  = X(1J)  + p(y+1)  + q(w+1) 


if  ||  x(p+1)  - X(U) 


is  less  than  some  prescribed  tol* 


erance  set  Yg(0)  = X^y  otherwise  increment  y by  1 
and  repeat  steps  (ii)  through  (vi). 

Of  course,  if  the  problem  is  almost  linear  then  only  one 
iteration  need  be  performed. 

The  entire  procedure  was  successfully  applied  to  several 
examples,  some  of  which  are  discussed  in  the  next  section. 

6.  NUMERICAL  EXAMPLES  AND  DISCUSSIONS  OF  RESULTS 

In  this  section  we  present  and  discuss  the  results  of 
three  examples  comparing  our  method  of  Section  5 with  Hind- 
marsh's  (1974)  version  of  Gear's  code  for  stiff  differential 
equations.  Both  the  Adams'  methods  that  we  use  to  integrate 
the  reduced  differential  equation  and  Gear's  stiffly  stable 
methods  are  contained  in  this  code,  and  the  user  sets  a 
parameter  to  select  the  appropriate  method.  Hindmarsh's 
code  and  the  IMSL  Bulirsch  and  Stoer  code  also  require  the 
user  to  select  an  estimate  for  the  relative  local 
discretization  error  and  an  initial  step  size  for  the  integ- 
ration. In  all  cases  we  selected  the  relative  error  toler- 


ance as  10 


-S 


This  is  perhaps  a bit  too  severe  for  our 


methods,  because  if  the  problem  is  not  very  stiff  the  reduced 

solution  will  be  calculated  with  more  accuracy  than  neces- 

-4 

sary.  The  initial  step  size  was  selected  as  10  for  Adams' 
methods,  l/10e  for  Gear's  methods,  and  1 for  Bulirsch  and 
Stoer's  method.  We  found  that  the  IMSL  code  was  extremely 
sensitive  to  the  choice  of  the  initial  step  size  and  the 
times  for  the  integration  varied  quite  dramatically  depending 
on  this  choice.  Our  choice  of  unity  seemed  near  optimal  for 
the  problems  that  we  considered. 

In  the  tables  that  follow  five  numerical  solutions  are 
compared.  The  solutions  labeled  "asymptotic"  were  calculated 
by  our  method  without  explicitly  recognizing  e and  using 
either  the  Hindmarsh  (Adams)  or  the  IMSL  codes;  those  labeled 
"Gear"  were  solved  by  Hindmarsh' s (Gear)  code;  and  those 
labeled  "reduced"  were  calculated  by  our  method  with  e expli- 
citly set  to  zero.  Additional  headings  in  the  tables  are  as 
follows : 

e is  the  maximum  difference  in  any  component,  times  106, 
between  a numerical  solution  and  the  exact  solution  at 
terminal  time  T.  For  our  asymptotic  or  limiting  solu- 
tions 

e = max  |Y  (T)  - y (T) | x 10&. 
l^i^m  u 1 

In  general,  for  the  examples  considered,  the  error  was 
fairly  constant  outside  of  the  initial  boundary  layer. 

d when  the  exact  solution  is  not  known,  we  have  tabulated 

6 

the  maximum  difference  in  any  component,  times  10  , 
between  solutions  obtained  by  our  method  and  those  by 
Gear's  code  at  terminal  time  T. 


NFE  For  our  asymptotic  or  limiting  solutions  this  denotes 
the  number  of  times  that  YQ  was  evaluated  during  the 
course  of  the  integration.  For  Gear's  solutions  it 


denotes  the  number  of  times  that  y was  evaluated. 

NJE  For  Gear's  solutions  this  denotes  the  number  of  times 


that  f was  evaluated  during  the  course  of  the  integ- 


Our  methods  evaluate  f each  time  Y„  is  eval- 

y o 


ration, 
uated. 

CP  Time  in  milli-seconds  to  integrate  the  problem,  exclu- 
ding input/output  and  supervisor  state  time.  Except 
where  noted  it  includes  the  time  necessary  to  calculate 
the  initial  conditions  YQ(0)  by  our  method.  In  most 
cases  the  times  are  averaged  over  several  runs.  All 
calculations  were  performed  on  an  IBM  360/67  at  the 
Rensselaer  Polytechnic  Institute. 

. CP  time  relative  to  the  fastest  execution  time, 
rel 

The  individual  examples  are  discussed  below. 

Example  1 : 

s-1)  2(1/ e-1)”]  r 1 1 

0 < t < T = 1. 

This  constant  coefficient  example  is  an  adaptation  of  one 
considered  by  Gear  (1971).  The  exact  solution  is 
-t  - -t/e 


CP 


"(1/e-l) 

2(1/ e-l)~ 

” 1 

y,  y(0)  = 

_-(l/e-l) 

-(2/e-l)_ 

_ 1 

y(t) 


4e 


-2e-t  +3e 


-3e 
-t/e 


i 


The  results  are  compared  in  Table  1 for  e = 10  , i=2,4,6,8. 

They  are  typical  of  the  results  of  subsequent  examples  in 
that  they  show  that  the  accuracy  of  our  method  increases  as 
the  stiffness  increases  without  an  increase  in  computational 
effort.  On  the  average,  our  asymptotic  and  reduced  solutions 
required  5 milli-seconds  to  calculate  the  initial  conditions 
for  the  outer  problem. 


TABLE  I 


\ 
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Summary  of  Results  for  Example  1 


e 

METHOD  10"2  10*4  10"6  10‘8 


Asymptotic 

e 

11100. 

111. 

1.81 

.0130 

(Adams ) 

* 

NFE 

34 

34 

34 

34 

CP 

135 

137 

138 

129 

CP  , 
rel 

1.57 

1.60 

1.62 

1.51 

Asymptotic 

e 

11100. 

110. 

1.10 

.00171 

(IMSL) 

NFE 

33 

33 

33 

33 

CP 

89.1 

89.8 

91.0 

89.8 

CPrel 

1.04 

1.05 

1.05 

1.05 

Reduced 

e 

.126 

.701 

.701 

.701 

(Adams) 

NFE 

34 

CP 

135 

C?rel 

1.58 

Reduced 

e 

.000309 

.000309 

.000309 

.000309 

(IMSL) 

NFE 

33 

CP 

85.5 

CP  . 
rel 

1.00 

Gear 

e 

3.51 

7.61 

8.85 

4.40 

NFE 

158 

183 

188 

195 

NJE 

17 

24 

25 

27 

CP 

332 

396 

407 

422 

CP  . 
rel 

3.89 

4.63 

4.75 

4.94 

TABLE  II 


Sumary  of  Results  for  Example  2 


e 

METHOD 

io'2 

io'4 

io'6 

io"8 

Asymptotic 

e 

300. 

1.98 

.640 

.666 

(Adams) 

NFE 

46 

46 

46 

46 

CP 

196 

192 

188 

196 

CP  . 
rel 

1.36 

1.33 

1.31 

1.37 

Asymptotic 

e 

301. 

2.65 

.0284 

.00223 

(IMSL) 

NFE 

49 

49 

49 

49 

V 

CP 

154 

150 

152 

148 

• i , * * •* 

CP  . 
rel 

1.07 

1.05 

1 .06 

1.03 

Reduced 

e 

214. 

2.81 

.688 

.667 

(Adams) 


NFE  46 
CP  187 
rp 

rel  1.30 


Reduced 

(IMSL) 

e 

NFE 

CP 

CP  . 
rel 

214. 

49 

144 

1.00 

2.14 

.0201 

.00175 

Gear 

e 

.541 

1.26 

3.93 

3.97 

NFE 

167 

191 

196 

203 

NJE 

20 

25 

26 

28 

CP 

341 

399 

406 

421 

CP  . 
rel 

2.38 

2.78 

2.83 

2.93 

TABLE  III 


Sumary  of  Results  for  Example  3 


METHOD 

-2 

10 

£ 

-4 

10 

-6 

10 

-8 

10 

Asymptotic 

d 

317. 

3.39 

.352 

.358 

(Adams) 

NFE 

30 

30 

30 

30 

CP 

170 

170 

170 

172 

CP  . 
rel 

1.66 

1.66 

1.66 

1.68 

Asymptotic 

d 

317. 

3.56 

.440 

.446 

(IMSL) 

NFE 

21 

21 

21 

21 

CP 

108 

107  107 

107 

•if  * ‘ 

CP  . 
rel 

1.05 

1.05  1.05 

1.05 

Reduced 

d 

1217. 

12.0 

.269 

.358 

(Adams) 

NFE 

30 

CP 

165 

CP  . 
rel 

1.61 

Reduced 

d 

1217 

12.1 

.411 

.445 

(IMSL) 

NFE 

21 

CP 

102 

CP  , 
rel 

1.00 

Gear 

NFE 

143 

144 

150 

158 

NJE 

17 

19 

21 

23 

CP 

343 

365 

380 

400 

CP  . 
rel 

3.35 

3.57 

3.72 

3.91 

l ___z n j 


TABLE  IV 


Time  to  integrate  from  t=0  to  t = T = 10c  for  Example  3 


e 

Asymptotic 

(Adams) 

Gear 

NFE 

CP 

NFE 

NJE 

CP 

d 

CP 

ratio 

CVJ 

1 

O 

12 

82.9 

119 

14 

316 

451. 

3.81 

.4 

10 

5 

52.1 

121 

15 

341 

4.81 

6.54 

o o 

i i 

co  a 

2 

37.1 

121 

15 

342 

.0361 

9.22 

2 

37.4 

121 

15 

341 

.0231 

9.12 

TABLE  V 


Time  to  integrate  from  t=0  to  t = T = 1 for  Example  3 
using  initial  conditions  for  the  outer  problem  ( results 
for  e = 0 are  the  reduced  solution) 


e 

Asymptotic/ 

Reduced 

(Adams) 

CP 

NFE  CP  rel 

Asymptotic/ 

Reduced 

(IMSL) 

CP 

NFE  CP  rel 

NFE 

Gear 

NJE  CP 

CP 

rel 

io“2 

30 

142 

1.87 

21 

79 

1 .04 

54 

6 

103 

1.36 

.4 

10 

30 

141 

1 .87 

21 

79 

1.04 

56 

5 

102 

1 .35 

• o 

10 

30 

141 

1 .87 

21 

• 79 

1 .04 

43 

8 

100 

1.32 

- 8 

10 

30 

143 

1 .89 

21 

78 

1 .04 

51 

10 

118 

1.55 

0 

30 

138 

1 .83 

21 

76 

1.00 

Example  2 : 

• = i f {yi+yZ)l1  ' I(yi+y2)2]  " 72  (y2"yi} 
e (yx-+y2) [i  - |(y!+y2)2]  + ~ 


0 S t < T = 2 


y(0) 


-2 

0 


This  nonlinear  problem  was  contrived  so  that  the  orthogonal 
matrix  E is  constant  and  the  exact  solution  is  known  as 
y 1 (t)  = (S  - n)/^  y2(t)  = (£  + n )/»/2 

with 


1 -2t/e.-l/2 
-(1  - 2e  ) 


/2 


1-1/C  \ -e 

Tin ' 


The  results  are  presented  in  Table  2 and  generally  parallel 
those  for  Example  1.  The  average  time  required  to  calculate 
the  initial  conditions  for  the  asymptotic  and  reduced  solu- 
tions was  24  milli-seconds 
Example  3: 


' (y|-y1y3)-ey' 

1 

9 = f(y,e)  - i 

2(y1y3-y^)+ey 

,y(0)  = 

0 

(y|-y1y3) 

1 

0<t<T  = 1 


This  example  arises  in  chemical  reactions  and  was  studied  by 
Vasil 'eva  (1976).  She  did  net  specify  the  eyj  terms  in  f 
nor  the  initial  conditions  and  they  were  selected  by  us 
rather  arbitrarily.  The  Jacobian  fy(y,0)  of  this  system  has 
rank  1 for  all  y =/  0 and  it  may  be  row-reduced  by  a constant 
orthogonal  matrix  E.  The  results  of  this  example  are  pre- 
sented in  Table  3.  The  average  time  required  to  calculate 
the  initial  conditions  for  the  asymptotic  and  reduced  solu- 
tions was  28  milli-seconds. 

Our  method  is  to  be  used  on  problems  where  the  boundary 
layer  solution  is  not  of  interest;  hence,  we  should  be  able 
to  calculate  the  initial  conditions  for  the  outer  problem 
faster  than  a stiff  differential  equation  solver  could  in- 
tegrate through  the  boundary  layer.  In  order  to  provide 
some  evidence  that  this  is  the  case  we  solved  Example  3 in 
the  interval  0 < t £ lOe  (the  approximate  boundary  layer 


region)  using  Gear's  methods  and  our  asymptotic  method  with 
the  Adams'  integrators.  The  results  are  presented  in  Table 
4 for  e = 10  i = 2, 4, 6,8.  The  CP  times  for  our  method 
includes  both  the  times  to  calculate  the  initial  conditions 
and  to  integrate  the  outer  problem  from  t = 0 to  10c.  To 
make  the  comparison  somewhat  more  fair  we  re-evaluated  E 
after  each  iteration,  even  though  it  is  constant  for  this 
example.  For  e ^ 10  we  see  that  our  method  can  calculate 
the  solution  at  the  edge  of  the  boundary  layer  region 
approximately  9 times  faster  than  Gear's  methods. 

A comparison  of  the  results  in  Tables  3 and  4 shows  that 
about  90%  of  the  time  required  to  integrate  Example  3 from 

t = 0 to  1 by  Gear's  code  is  devoted  to  the  boundary  layer 

-4 

region  for  e ^ 10  . This  suggests  the  possibility,  of  using 

our  method  to  calculate  the  initial  conditions  for  the  outer 
problem  and  then  using  a stiff  method  to  integrate  the  ori- 
ginal differential  equation.  This  test  was  performed  on 
Example  3,  and  the  results  are  reported  in  Table  5.  All 
methods  use  the  same  initial  conditions,  i.e.,  those  gener- 
ated by  our  method.  The  CP  times  required  to  calculate 
these  conditions  are  not  included  in  Table  5.  The  differ- 
ence between  any  two  computed  solutions  is  less  than 
-4 

3 x 10  . While  the  results  are  far  from  conclusive,  they 

do  show  the  extra  computational  effort  that  is  required  by 
Gear's  method  for  very  stiff  problems. 

The  state  of  the  art  of  numerical  methods  for  stiff 
initial  value  problems  for  ordinary  differential  equations 
is  very  well  developed  (cf.  Enright  at  al  (1975))  and  a 
variety  of  good  techniques  exist.  Nevertheless,  there  are 
many  problems,  particularly  in  chemical  reactions,  where 
asymptotic  methods  should  be  useful.  They  may  be  used  to 
calculate  accurate  solutions  of  very  stiff  problems,  to 


furnish  initial  conditions  for  standard  stiff  integration 
routines,  and/or  as  an  analytical  tool  to  provide  qualitative 
information  about  the  solutions  of  stiff  problems.  In  future 
papers  we  hope  to  extend  our  calculations  to  initial  value 
problems  where  E depends  on  y and  to  consider  boundary  value 
problems. 
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